function [diff, eq] = find_stationary_dist(eq, param,glob,options)
    
    Q           = eq.Q;

    % Distribution of incumbent firms
    L               = eq.L;  
    
    for itL = (1:options.itermaxL)
        Lnew    = Q'*L;  
        dL      = norm(Lnew-L)/norm(L);
        if (dL<options.tolL),break,end

        L       = Lnew;
    end
    
    L = param.M*L;

    % L is the distribution over all states
    % Lp is the distribution over customer bases

    Lp = nan(size(glob.bgridf));
    for i = 1:length(glob.bgridf)
        Lp(i) = sum(L(find(abs(glob.sf(:,1) - glob.bgridf(i)) < 1e-6)));
    end

    eq.L  = L;
    eq.Lp = Lp;
    
    q = eq.yf/param.C;
    diff = (sum(eq.L.*upsilonp(q,param).*q))^(-1) - param.D; 

    eq.Dout = (sum(eq.L.*upsilonp(q,param).*q))^(-1);
    
end